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Abstract 


The present research is concerned with 2-D simulations of turbulent flows using Spalart- 
Allmaras turbulence model on various geometries like flat plate, NACA0012 airfoil, Mc- 
Donnell Douglas SOP-SON multi-element airfoil. The incompressible Reynolds Averaged 
Navier Stokes(RANS) equations, in conjunction with the Spalart-Allmaras turbulence 
model for turbulence closure, in primitive variables, have been solved using stabilized 
finite element formulations. Typical finite element mesh consists of structured mesh 
close to the body and an unstructured mesh, generated via, Delaunay’s triangulation, 
in the rest of the domain. This type of grid has the ability of handling fairly complex 
geometries while still providing the desired resolution close to the body, to effectively 
capture the boundary layer and shear layer, especially in the context of unsteady flows. 
In addition, the structured grid around the body allows for efficient implementation of 
the turbulence model. 

The results have been presented for turbulent computations. Flat plate results have 
been found to be in good agreement with the theoretical results with successful control 
over the trip locations. The flow phenomenon over NACA0012 airfoil has also been 
captured successfully. In the context of multi-element airfoil, the computed results have 
been presented at Re — 9 x 10 6 , for angle of attack varying from 8° to 27°. It has been 
found that the velocity profiles match quite well with the experimental results below the 
slat wake but gets diffused towards the downstream of the flow even at lower angle of 
attack. C p distribution also matches quite well on the lower surface but shows a slight 
deviation on the upper surface at the higher angle of attack. Since the turbulence model 
is not able to capture the stall behaviour at higher angle of attack, the lift coefficient 
has been found over predicted in computations. The computations of slat flow field 
represents a key roadblock for successful prediction of multi-element flows. 
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Chapter 1 


Introduction 


High-lift aerodynamics has received considerable attention in last few years. A number 
of conferences and workshops have been held on high-lift problems to search for better 
computational and experimental tools. Also significant research works are in progress to 
resolve the difficulties, while the prediction of high-lift flow-fields consisting of multiple 
airfoil elements still remains a challenge to the computational fluid dynamics (CFD) 
community. Multi-element airfoil has wide variety of flow physics which requires accu- 
rate modeling. Each airfoil element has a stagnation point and boundary layers that 
developed along a curved surface in pressure gradients and may pass through various 
transition of laminar to turbulent flow. Additionally, the wake of each element eventually 
merges with the boundary layer, developed from subsequent elements. Hence keeping 
the slat wake distinct from the main element boundary layer for the greatest distance 
on the main element is one of the key ingredients of high lift multiple element design. 

The factors affecting the design of multi-element airfoil namely, pressure distribu- 
tion over surfaces and the flow around the slat, is also quite complex. The upper surface 
of the slat is highly curved with large flow acceleration, whereas the flow at the slat trail- 
ing edge can be laminar, transitional, or turbulent. On the other hand, the slat lower 
surface is highly concave and a result of that the flow in the lower side is accelerated 
through the slat-main element gap around a very sharp bend. The flow often separates 
due to the slat geometry and reattached near the slat trailing edge. The free shear layer 
developed at the recirculation boundary is most likely turbulent. The wake, developed 
at the trailing edge of the slat, is created from the laminar, transitional, or turbulent 



1.1 Problem description 


2 


upper surface flow and the turbulent lower surface flow. Additionally, the flow around 
the lower surface may be unsteady. 

Therefore, the predictive capability for high lift multi element airfoil configurations 
is dependent on the accurate modeling of (i) boundary layer transition on a curved surface 
in a strong pressure gradients (ii) turbulent boundary layer development on a curved 
surface in a strong pressure gradient (iii) transitional and turbulent wake development 
(iv) separating and reattaching flow and (v) wind tunnel effects. 

To achieve the goal of a high-lift system, generating as much lift as possible, the 
flow should remain attached to the airfoil surface. This can be achieved effectively 
by using multiple elements instead of using any external devices such as wall suction. 
Multi-element geometry provides a greatest advantage of getting high-lift through invis- 
cid pressure distribution over surfaces by reducing the pressure rise over each element. 
However, the presence of multiple elements complicates the analysis procedures because 
of important and often complex interactions between the individual elements. While in- 
viscid analysis can be done using panel methods or unstructured grid Euler solvers. But 
the accurate prediction of the flows about these configurations requires viscous analysis 
which can have large influences on the pressure distributions. As a result, these have 
obvious effects such as separation on the surfaces, wake interactions between forward 
and aft elements as well as flow reversal off the surface, can contribute in determining 
the overall performance of the high-lift system significantly. 

The simulation of real flows, viscous in nature, has become a complicated task and 
hence the solution of Navier Stokes equation is essential to capture the effects of viscosity 
such as skin friction drag and heat transfer coefficient. Furthermore, in the real world 
problems involving high Reynolds number as those in aerospace applications, the flow is 
turbulent in nature. 


1.1 Problem description 


Multi-element airfoil has important applications in high-lift systems, typically for present 
transport airplanes, often consisting of a leading edge slat and a trailing edge flap system. 
The analysis brings forth unexpected complexities in flow structures. The objective 
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of the present thesis is to develop a tool for analysis of such flows and their physical 
understanding. 

The present study is concerned with the high Reynolds number flows obtained from 
two-dimensional numerical computations. The unsteady incompressible Reynolds Av- 
erage Navier Stokes (RANS) equations, in conjunction with Spalart-Allmaras model 
for turbulence closure in primitive variable formulation have been solved on hybrid 
mesh with triangular elements using a well proved stabilized finite element method with 
the SUPG(streamline-upwind/Petrov-Galerkin) and PSPG(pressure-stabilizing/Petrov- 
Galerkin) stabilization terms. The resulting nonlinear equation system has been solved 
using Generalized Minimal RESidual (GMRES) technique in conjunction with diagonal 
pre-conditioners . 

Results have been compared with other published results in the open literature. 
The effect of spatial resolution has also been investigated. The flow problem considered 
in the present thesis are turbulent computations over (a) flat plate, (b) NACA 0012 
airfoil (Re = 10 6 ), and (c) 30P-30N multi-element airfoil (Re = 9 x 10 6 ). 


1.2 Literature survey 


The literature pertaining to flow past multi-element airfoil is reviewed below. The review 
of earlier investigations related to the different aspects of the flows for high Reynolds 
number has also been presented. A brief review on the modeling aspects of turbulent 
flow has also been undertaken to examine the performance of various turbulence models 
for complex geometries. 

Some important insight into the physical processes that govern the high-lift aerody- 
namics were summarized by Smith in his landmark paper [1]. In particular he described 
several effects that contribute to the improvement of high lift characteristic of multiple 
elements over single element. 

Turbulence modeling is an important issue in solving flows past complex geome- 
tries. The various length scales must be resolved properly. Tulapurkara [2] gave an 
extensive review on the development and application of turbulence models to aerospace 
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applications. Spalart-Allmaras [3] proposed a one-equation turbulence model. This is 
one of the most widely used model for aerodynamic flows in recent times. They came 
out with some advantages over Baldwin-Barth’s model. Thereafter, Spalart-Shur [4, 5] 
came out with a proposal for sensitization of the model to rotation and curvature effects. 

Anderson et. al. [6] utilized a two-dimensional unstructured Navier-Stokes code for 
computing the flow around multi-element configurations using both Baldwin-Barth and 
Spalart-Allmaras turbulence models. Comparisons had been shown for a landing config- 
uration with an advanced-technology flap. They also showed the comparisons between 
experimental & computed results and the grid convergence studies to assess inaccuracies 
caused by inadequate grid resolution. They concluded from the grid convergence study 
that sufficient grid refinement is required to obtain suitable levels of grid convergence for 
velocity profiles at high angle of attacks, while adequately small spacing of grid points to 
be maintained for obtaining accurate resolution of the wall boundary layers on upstream 
elements. 

Godin et. al. [7] reported detailed comparison of the Spalart-Allmaras and Menter 
turbulence models in the context of two-dimensional high-lift aerodynamic flows. Their 
results showed that the Menter model to be more accurate in separated flow regions, 
whereas the Spalart-Allmaras model is more accurate in attached flows and wakes, in- 
cluding merging boundary layers and wakes. The Spalart-Allmaras model is somewhat 
more robust, though for several cases the computational costs were about equal. 

Rango and Zingg [8] investigated the flow over multi-element airfoil configurations 
using higher-order methods on multi-block grids, in conjunction with Spalart-Allmaras 
turbulence model. They performed grid convergence study to evaluate the accuracy of 
the higher-order algorithm on multi-block grids and showed that the globally higher- 
order approach can provide accuracy equivalent to a second-order algorithm on a grid 
with several times fewer nodes, thus greatly reducing the cost of computing flows over 
multi-element airfoils. 

Rumsey and Gatski [9] applied the Spalart-Allmaras model with curvature correc- 
tions to simulate the flow over multi-element airfoil. Their investigations showed that 
implementation of curvature term does not have much impact on flow field. 

Rumsey et. al. [10] demonstrated the need to include the transition location in 
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addition to study of various turbulence models. They studied the McDonnell Douglas 
three element high lift configuration tested in NASA Langley LTPT(Low turbulence 
Pressure Tunnel) wind tunnel. They derived a conclusion from their study that transi- 
tion location is crucial to the accurate computation of boundary layer velocity profiles. 
Otherwise poor boundary-layer profile predictions lead to poor predictions of trailing- 
edge near-wake profiles and, therefore, poor predictions in the far-wake region. They also 
demonstrated their efforts in understanding the flow over slat by specifying transition on 
the slat resulted either in slat separation or in only marginal improvement in slat-wake 
profiles compared to the experimental results. 

Kim and Nakahashi [11] solved RANS on the unstructured grid for computing 
flows over multi-element airfoil configuration. They compared their results with the 
experimental data to assess the accuracy of the computations. They also discussed the 
effect of complex flow phenomenon around high-lift system, like the interaction between 
upstream element wakes and boundary layers on downstream elements, on the solution 
accuracy using the grid refinement method. 

Valarezo and Mavriplis [12] presented the application of compressible RANS method 
to the calculation of flows about a transport-type multi-element airfoil. They used 
unstructured-mesh method utilizing multi grid techniques for computational efficiency in 
the context of different turbulence models. They discussed the usability of such method 
to capture performance increments due to rigging changes and Reynolds number effects 
and suggested to consider grid adaptation for extensions to three-dimensional flows. 

Joseph H. Morrison [13] performed the calculations with Reynolds averaged Navier- 
Stokes code ISAAC to study the effect of numerical accuracy on the advection terms of 
the turbulence equations and the performance of several turbulence models. Also the 
effect of transition location was investigated. His results estimated too large wake defect 
in all the models and finally suggested to go for improved transition model. 

Jones et. al. [14] discusses the multi zonal grid technique for analyzing the multi- 
element airfoil. Grid resolution is an important aspect to capture the shear layer of each 
element and for that they conducted grid refinement study and found dramatic effect on 
the downstream flow field velocity profiles. 

Rogers et. al. [15] reported the comparison of various turbulence models in com- 
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puting multi-element airfoil flows with the experimental results. Their results showed 
very little difference in C p distribution in the context of various models. They also re- 
ported that Spalart-Allmaras(SA) model appeared to have the most amount of mixing 
in the wakes, while Shear stress transport(SST) model consistently showed the greatest 
amount of velocity defect in the wakes. According to their results, SA model was the 
only one to show a stall behaviour below an angle of attack of 25 degrees, this may 
be due its treatment in transition. The most significant thing observed in their study 
is that the difference between the computational and the experimental data in the slat 
wake at lower angle of attack. The grid resolution study showed that this difference did 
not decrease with increased grid resolution. They finally concluded that all the models 
are similar for most part of the predictions but given the complexity of the multi-element 
airfoil geometry and the surrounding flow field, further development and tuning of tur- 
bulence models were suggested. Nakayama et. al. [16] did experimental investigation of 
flow field about multi-element airfoil to measure the mean flow and turbulence quantities 
using pressure and hot-wire probes. 

Experiments have also been performed in NASA LTPT on McDonnell Douglas 
three element landing configuration, which was used as a test cases in a CFD challenge 
workshop discussed in [17] held at NASA. The conclusions drawn on the workshop 
were (i) RANS showed less variability than did potential/Euler solvers coupled with 
boundary-layer solution techniques, (ii) Coupled-methods drag prediction agreed more 
closely with experiment than the RANS methods. Lift was more accurately predicted 
than drag in both methods, (iii) Both the methods did reasonably well in predicting lift 
and drag due to changes in Reynolds number, (iv) Pressure and skin friction were also 
in good agreement with the experimental results, and finally it was suggested to have 
accurate transition locations since it appeared to have a major influence on overall airfoil 
performance. 


1.3 Thesis organization 

Reynolds Averaged Navier-Stokes equation have been presented in chapter2, followed 
by an introduction to turbulence modeling. The Spalart-Allmaras turbulence model has 
also been described in the same chapter. The details of numerical techniques, governing 
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equations and formulations with the different boundary conditions including the grid 
system used here have been presented in Chapter3. Chapter4 describes the results of 
the present research and related discussions. Chapter5 summarizes major conclusions of 
the present work. The scope for future work has been enlisted in Chapter6. 



Chapter 2 


An Overview of Turbulence 
Modeling 


Newtonian fluid flow is now understood to be fully represented by the Navier-Stokes 
equations. Numerical solution of Navier-Stokes equations is thus a successful modeling 
tool for laminar and turbulent flows. If turbulent flow is computed using the same grid 
and time as in laminar flow, the solutions will not reveal the existence of the small scale 
eddies. Hence, the transport of turbulent kinetic energy from the larger to the smaller 
eddies will not be addressed. Nevertheless, the Navier-Stokes equations can be solved on 
a fine enough grid with an exceptionally accurate discretization method so that both fine 
scale and large scale aspects of turbulence can be calculated. This approach is termed as 
the Direct Numerical Simulation(DNS ) of turbulence. DNS has been a very successful 
tool over the past ten years for the study of turbulent flow physics. It has however 
a severe limitation. In order to resolve all scales of motion, one requires a very large 
amount of computing resources. 

In Large-Eddy Simulation^ LES) the small scales of turbulence, which are assumed 
to be more universal in nature, are modeled; the large, energy-containing scales that 
are more flow specific are explicitly calculated. The models used for the small scales 
are called sub grid-scale(SGS) models. Unlike other types of turbulence models, the 
resultant equations for LES describe a fully time-dependent flow; the modeling only 
blurs the turbulent structures so that the small scales do not need to be calculated. Like 
the DNS, the LES also requires enormous computing resources for solving practical flow 
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problem. A potential exists for the practical use of LES & DNS in the near future, 
although the high cost will restrict its use to cases where lower order models are not 
expected to give satisfactory results. 


2.1 Statistical Turbulence Models 

Presently, statistical turbulence models are the most widely used turbulence-modeling 
schemes. Flow quantities are decomposed into mean and fluctuating parts and then 
substituted into the equations of motion. These equations are averaged to produce a 
set of equations for the mean motion. Osborne Reynolds [24] pioneered this approach 
by temporally averaging the equations for incompressible flow; the resultant equations 
are known as the Reynolds-averaged Navier-Stokes(RANS) equations. These equations 
involve the mean flow quantities as well as correlations of the fluctuating quantities. 
The correlations appear in the equations of mean motion in the same way as the viscous 
stress terms appear; hence, these correlations are known as the Reynolds stresses. The 
various classes of models differ in how these correlations are approximated. 

Reynolds stress transport (RST) models involve transport equations for each of the 
six independent Reynolds stress components. This class of models is the most complex of 
the statistical turbulence models, and the use of these models for engineering applications 
is not yet commonplace. However, because the Reynolds stresses can independently 
respond to various flow conditions, this class of models can potentially be applied to a 
large variety of flows. This potential generality motivates much of the current research 
on these models. 

Eddy-viscosity models include a number of classes of models, all of which approx- 
imate the effect of the turbulence on the mean motion by modifying the coefficient of 
viscosity. The effective viscosity coefficient that issued in the computation of the flow 
field is the sum of the molecular viscosity(/i) and turbulent viscosity (fa). The different 
classes of eddy- viscosity models are distinguished by the number of additional differential 
equations that are solved to determine fa. Dimensional analysis suggests that fa is the 
product of the density, a velocity scale, and a length scale. The local mean density is 
almost always used, leaving the velocity and length scale still to be determined. Two- 
equation models solve differential equations to determine these two scales. One-equation 


2.2 Integral Methods 


10 


models solve a differential equation for the velocity scale and use algebraic relations to 
determine the length scale. Zero-equation models use algebraic relations to determine 
both the velocity and length scales. More detail on all of these model classes are given 
in books. 


2.2 Integral Methods 


In integral methods, an ordinary differential equation is solved for the momentum thickness (0) 
in terms of the skin-friction coefficient (Cf), the displacement thickness(<5*), the boundary- 
layer edge velocity, and body curvature. For flows with heat transfer, another equation is 
required that involves the Stanton number (St), the enthalpy thickness, the free-stream 
and wall temperatures, and the body geometry. Approximate relationships between 
the variables are substituted into the equations. The equations are integrated in the 
downstream direction. These methods are computationally quite efficient, but are ac- 
curate only for those flows in which the assumed relationships are appropriate. Of the 
20 integral methods that competed in the 1968 Stanford Conference ( [25], [26]), only 
four received an evaluation of ’’good” after they were tested for 16 different turbulent 
flows. In spite of their limited generality, integral methods provide good skin-friction 
and heat-transfer predictions for well-studied flows of engineering interest. This method 
is discussed in more details in book [38]. 

2.3 Brief description of turbulence 
equations (RAN S ) 

The unsteady Navier-Stokes equations are generally considered to govern the turbulent 
flows in the continuum regime. Turbulent flow, shows irregular fluctuation in the flow 
variables, are superimposed on the mean flow. Therefore, the mathematical descrip- 
tion of turbulent flow, consists of decomposing the flow variables like velocity, pressure, 
density, temperature etc. into mean and fluctuating components (called as Reynolds de- 
composition) are defined as follows. 
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Now by time averaging the N-S Equations, and taking into account the basic rules for 
fluctuating components, we obtain RANS Equations: 



(2.3.5) 


l ( “ f) + ^ {Vi) = -^ + s: (T * ‘ (2 - 3 - 6) 

For the present cases, in the limit M — > 0 the governing equations are incompressible 
in nature i.e. density p is constant and the energy equation is decoupled. Hence, the 
continuity and momentum equations become 
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(2.3.7) 

(2.3.8) 


(2.3.9) 


Thus the mean momentum equation is complicated by new term involving the turbulent 
inertia tensor This new term may not be negligible in a turbulent flow, in general, 
and is the source of analytic difficulties. The components of u' { u are related not only 
to fluid physical properties but also to local flow conditions (velocity, geometry, surface 
roughness, and upstream history). If we observe equation 2.3.8 then, a turbulent inertia 
term behaves as if the total stress on the system were composed of the Newtonian 
viscous stresses plus an additional or apparent turbulent stress tensor which arises from 
the fluctuating terms. These turbulence-correlation terms must be modeled with semi 
empirical assumption, which leads to idea of turbulence modeling. With the above 
background information, the turbulence models are categorized. 
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2.4 Classification of Turbulence models 


The turbulence models can be classified depending on the number of differential equations 
solved in addition to the mean flow equations: 

1. Zero equation models - Algebraic models. 

2. One equation models. 

3. Two equation models. 

4. Stress equation models. 

Models in class 1-3 are based on Boussinesq eddy viscosity model, while those in 4 
obtain the Reynolds stress solving differential equation for u'v' directly. Zero equation 
models, which uses only the partial differential equation for the mean flow field and 
no transport equations for turbulence quantities, is also called ’’mean field” closure. 
The class 2-4 are called ’’transport equation” closures, as they takes into account the 
transport phenomena involved. 


2.5 Desirable features of Turbulence model 


The turbulence models provide the equations necessary to make the Reynolds Averaged 
Navier Stokes equations a closed set. The turbulence model should have the following 
features : 

• predict accurately a wide range of flows; 

• be mathematically simple and involve small number of assumptions and model 
constants; 

• be computationally stable and economical in computational time. 

• be dimensionally correct. 
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2.6 Spalart-Allmaras Turbulence Model 


This model was proposed by Spalart and Allmaras [3]. The proposition of this model 
was prompted by Baldwin and Barth’s [27] work, and by the belief that generating a 
One-equation model as a simplified version of the k — e model is not optimal. A one- 
equation model is simple enough that it can be generated ’’from scratch”, which may 
lead to better performance and certainly gives fuller control over its mechanics. A case 
in point is the Baldwin-Barth diffusion term, which constrained by the k — e ancestry 
and the further assumptions made. This model allows a ’’semi-local” near- wall term, 
also has the same properties as that of Baldwin-Barth in terms of compatibility with 
unstructured grids and benign near-wall behavior, and is more accurate, especially away 
from wall, as well as slightly more robust. 

The model has four nested versions from the simplest, applicable only to free shear 
flows, to the most complete, applicable to viscous flows past solid bodies and with 
laminar regions. As each additional physical effect is considered, new terms or factors 
are introduced. They are identified by a common letter subscript in the constants and 
functions involved. The new terms are passive in all lower versions of the model, so that 
the description proceeds in order. 


Free shear flows 

The central quantity is the eddy viscosity v t \ the Reynolds stresses are given by the 
constitutive relation —UiU] = 2u t Sij. For generating transport equation for v t at high 
Re, the molecular viscosity is not incorporated here. The left hand side of the equation 
is the Lagrangian or material derivative of u t \ Dv t /Dt = dv t f dt + Uidu t /dxi and on the 
right hand side the production and diffusion terms are provided. 

For the production term, the deformation tensor dUi/dxj presents itself. Since v t is 
a scalar, the scalar measure of the tensor is needed, denoted by S, defined as magnitude 
of the vorticity(| u> |), and then Su t has the desired dimension. The production term, 
and in fact the restriction of the model to homogeneous turbulence, is 

= c bi Su t (2.6. TO) 

The subscript b stands for ’’basic”. The eddy viscosity is stationary in isotropic turbu- 
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lence {i.e.Dv t / Dt = 0, since 5 = 0). The calibration of the model on inhomogeneous 
flows yields values of c b i between 0.13 and 0.14. 

The diffusion terms naturally focuses on spatial derivatives of v t . Incorporating 
the non-conservative diffusion term, involving first derivatives of v t , the ” basic” model 
defined as 

= c bl Su t + I [V. (v t Vv t ) + c 62 (Vi/ t ) 2 ] (2.6.11) 

where the constants are defines as a' = 2/3, c 61 = 0.1355, c b2 = 0.622. 

Near-wall region, high Reynolds number 


In a boundary layer the blocking effect of a wall is felt at a distance through 
the pressure term, which acts as the main destruction term for the Reynolds shear 
stress. This suggests a destruction term in the transport equation for the eddy viscosity. 
Dimensional analysis leads to a combination — Cwi(vt/d) 2 as a starting point, where d 
is the distance to the wall. The subscript w stands for a ’’wall”. This term is passive 
in free shear flows (d 5, so that the new term is much smaller than the diffusion 
term) and therefore does not interfere with the calibration upto this point, the idea 
of near- wall, but not viscous, ’’blocking’ term is also defined in other literatures. It is 
related to algebraic models, which takes the smaller of two eddy viscosities. In these 
models the outer eddy viscosity scales with the boundary-layer thickness, and the inner 
eddy viscosity is given by the mixing length, l oc d. In a classical log layer it is found 
that S = u T j (k d) and v t = u T nd follow these relationships. Equilibrium between the 
production and diffusion terms (all positive) and the destruction term is possible provided 
Cwi = c b i/t c 2 + (1 + c b2 ) /o'- 


Tests show that the model, when equipped the destruction term, can produce an 
accurate log layer in a U + (y + ) plot(this requires an adequate treatment of the viscous 
region, which is described later on). On the other hand it produces too low a skin- 
friction coefficient in a flat-plate boundary layer. This shows that the destruction term 
as formulated above decays too slowly in the outer region of the boundary layer. To 
address this deficiency and allow a new calibration, the term is multiplied by a non- 
dimensional function /„,, which equals 1 in log layer. With these modifications the model 
becomes 


= c bl Su t + ~ [V. iy&vt) + c b 2 (V^t) 2 


Cwlfw 



( 2 . 6 . 12 ) 
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The choice of an argument for f w was inspired by algebraic models, in which the 
mixing length plays a major role near the wall. This length can be defined by l = <J(v t /S) 
and the square of l/nd is considered for convenience: 


r = 




fw(r) = g 

with the constants as k 


1 + c! 


'wZ 


U 6 + 4s J 

0.41, Cyj2 := 0.3, Cyj 3 


Stfd 2 

/e 

, g = r + c w2 (r 6 -rj 


2.0 


(2.6.13) 

(2.6.14) 


Near-wall region, finite Reynolds number 


In buffer layer and viscous sublayer, additional notation is needed. Besides the wall 
unit, y + and so on, P is introduced which is equal to v t except in the viscous region. Also 
X follows the relationships^ = «y + from wall to the log layer. The transported quantity 
P behaves linearly near the wall. This is beneficial for numerical solutions: P is actually 
easier to resolve than U itself, in contrast with e, for instance. Therefore, the model will 
not require a finer grid than an algebraic model would. To arrive at this behavior, the 
classical law of the wall is considered and devised near- wall ” damping functions” which 
are distinct from the f w near-wall non-viscous destruction term. 


The eddy viscosity v t equals tzyu T in the log layer, but not in the buffer layer and 
viscous sublayer. P is defined in such a way that it equals Kyu T all the way to the wall. 
This leads to 3 

■* = »/. i. /»i = 3T7T (2.6.15) 

The subscript v stands for ’’viscous” and Cy i contains a value of 7.1. 

To modify the production term, S is replaced with S, given by 

S,S + ^U, /* = (2.6.16) 

The function / u2 is constructed, just like f v i, so that S maintains its log-layer behavior(5 = 
u T / («;y)) all the way to the wall. So, the other quantities involved in the ’’inviscid” model 
are redefined in the terms of P instead of u t , for instance r = P/ (Sk 2 cP'). Finally vis- 
cous diffusion term is added, consistent with a Dirichlet boundary condition at the wall, 
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v = 0. And the transport equation has become 

— = c bl Su + - [V. ((i/ + P) VP) + c 62 (VP) 2 ] - c wl f w [-J (2.6.17) 

This equation now yields equilibrium [Dv/Dt = 0) all the way d = 0 in a classical law- 
of-the-wall situation. 

Laminar region and trip term 

The final set of terms provide control over the laminar regions of the shear layers, 
which has two aspects: keeping the flow laminar where desired, and obtaining transition 
where desired. 

The turbulence index, i t = ( where u T = \J[y \ oj |)), determines the tran- 

sition near a wall. This index attains an approximate value of 0 & 1 in a laminar and 
turbulent region, respectively. In laminar region P is less than about v and the produc- 
tion term is multiplied by a factor (1 — fy), where 

ft2 = Ct 3 exp (-Ct 4X 2 ) (2.6.18) 

The subscript t stands for ’’trip” and the constants are having the value of c* 3 = 1.2, ct 4 = 
0.5. In all case c t3 must be larger than 1. Particularly for boundary layer calculation the 
f t 2 term can be left out by assigning = 0. 

To initiate transition near the specified trip points in a smooth manner, a source 
term needs to be added that will be non zero in a small domain of influence and this 
domain should not extend outside the boundary layer. Not wanting to find this edge, 
nor violating the invariance principles, two quantities (A[/, Lo t ) are invoked, where A U is 
the norm of the difference between the velocity at the trip(i.e. usually zero since the 
wall is not moving) and that at the considered field point. u> t is the vorticity at the wall 
at the trip point. Also d t , the distance from the field point to the closest trip point or 
line, is introduced. Dimensional analysis points to AZ7 2 as a starting point for the source 
term, and the final equation becomes 

^ = c bl [1 - f t2 ] 5P+^ [V. ((u + P) VP) + c b2 (VP) 2 ] - \c wl f w - ^/ t2 ] j^] 2 +/ tl AE/ 2 

(2.6.19) 
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with 

fa = ctigti exp ^-c t2 [d 2 + g 2 d 2 ] (2.6.20) 

and g t = min (0.1, U/u> t Ax t ), where x t is the grid spacing along the wall at the trip, 
whereas the constants Cn,Ct2 hold the value 1.0 & 2.0, respectively. The Gaussian in / tl 
confines the domain of influence of the trip terms as needed; it is roughly a semi-ellipse. 

In the present code, the trip terms have not been employed except flat plate 
calculations. So, the resultant equation effectively solved, is equation 2.6.17. 

Curvature effects 

The implementation of curvature terms to the one-equation Spalarat-Allmaras 
model lead to a fairly simple modification of the original model. The only difference 
between the modified model, accounting for the rotation-curvature effects and the SA 
model is that in the former the production term in the eddy viscosity transport equation, 
c b iSu, is multiplied by a rotation function / rl , defined as follows: 

2 r * 

fri ( r *, f) = (1 + Cri) — — [l - Cr Z tan~ x (Cr 2 f)] - Cn (2.6.21) 

Assuming that all of the variables and their derivatives are defined with respect to the 
reference frame of the calculation, which is rotating at a rate ft, the non dimensional 
quantities r*and f are given by the following formulas: 

no.. I 

r* = S/u, f = 2w it S 3 - l [- 5 ^ + (e lm „S J -„ + e jm „S j „)Sl m J/D 4 (2.6.22) 

Here, = 0.5 (f* 4 + £j) , = 0.5 [(^ - gf) + 2 £mji fl m ] 

S 2 = 2 SijSij, u? = 2 UijUij, D 2 = 0.5 (5 2 + u 2 ) 

and are the components of the Lagrangian derivative of the strain tensor. 
These axe all second derivatives. The Einstein summation convention is used with the 
model constants Cri = 1.0, Cr 2 = 12.0, = 1.0, respectively. 

In the present study, there is no rotation effect, i.e. Q = 0.0, so the 
curvature factor also gets simplified. It has also been observed that the implementation 
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of the curvature term leads to a slow convergence rate without having any appreciable 
change in the solution, so the curvature terms have been no more considered in the rest 
of the computations. 


2.6.0. 1 Advantages and Disadvantages of Spalart-Allmaras Model 

Some of the advantages of Spalart-Allmars model are: 

1. It is applicable to free shear flows and also to viscous flows past solid bodies with 
laminar regions. 

2. It is more accurate in attached flows and wakes, including merging boundary layers 
and wakes. 

3. It is also very much useful for computations in unstructured grid. 

4. In calculation of eddy viscosity, the upstream value are taken into account which 
improves the solution. 

Some disadvantages are: 

1. It is not a simple model for implementation, because it involves an extra differential 
equations and takes large computational time. 

2. It can only predict the velocity profiles u and v and the turbulent shear (— u'v') 
but cannot give the turbulent energy or any of the components ( u ', v', w') rms . 


Chapter 3 


Numerical Technique and Procedure 


The 2-D flow simulations over multi-element airfoil are carried out using incompressible 
Navier-Stokes equations solver, by stabilized finite element method. The finite element 
method starts with a piecewise approximation to the dependent variables. Various meth- 
ods of this class exist, all requiring an integral representation of the partial differential 
equation to be constructed. The classical finite element method for structural mechanics 
are based on variational principles. But for many engineering problems, particularly in 
fluid flow, more general approaches, such as the method of weighted residuals are used. 
The primary advantage of FEM is its ability to handle a complex geometry. The major 
activities involved in finite element analysis broadly classified into three classes: 

1. Pre-processing 

2. Processing 

3. Post-processing 

Pre-processing refers to all the tasks to be performed before starting the flow 
solution using N-S solver. The major pre-processing activities are smoothing the airfoil, 
generating the mesh and applying the required boundary conditions. The other impor- 
tant activity is to automate the grid generation procedure so that different grids can be 
generated for different configurations like airfoil and parafoil at different angles of attack. 
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Smoothing the Airfoil: The airfoil geometry which we obtained with contains 
very few points. These many points are not sufficient to generate a fine mesh which can 
resolve the flow field properly. So, to increase the number of points on the airfoil surface, 
we used a NAG routine E02BAF and E02BBF which fits cubic splines. 

Mesh generation: This is one of the very important task in CFD because with a 
poor mesh one can never expect good results even with higher order accurate method. 
The mesh generation described in the present study, follows Vinod’s work [28]. Mesh is 
basically set of points, with certain connectivity which one obtains by discretizing the 
domain of flow. Because of this spatial discretization the continuum of space is replaced 
by a finite number of points where numerical value of variables has to be determined. 
For the simulations the space domain has been discretized using hybrid grid i.e. the 
region close to the body consists of structured mesh while the rest is filled with unstruc- 
tured mesh. Unstructured meshes are very useful for the complex geometries at the same 
time they take into account the complexity of the flow. Structured meshes have logical 
indexing scheme for the nodes while, the unstructured meshes have a irregular connec- 
tivity. In the structured mesh generation the points are distributed in any ratio uniform 
or geometric progression or exponentially varying. In our cases, structured meshes are 
constructed using lines normal to the surface of the body. Triangular elements have been 
considered which are supposed to be ideal for the complex geometries. Filling up the 
whole domain with structured mesh would be very costly and controlling the density of 
points would be difficult. So to capture the shear layer coming out of slat we generated a 
structured mesh just ahead the slat. Further refinement is done by two level refinement 
of the main element mesh. Here the 2nd level of refinement is adjusted according to slat 
shear layer. Over the flap we inserted a inclined fictitious boundary following the flap 
deflection. 

For the unstructured mesh generation, among various algorithms available for tri- 
angulation, Delaunay approach has proved to be very efficient. The most important 
property of a Delaunay mesh is that for each element in the mesh, its circumcircle which 
encompasses all nodes defining that element contains no other nodes of the mesh. This 
produces elements with very optimal aspect ratio. The algorithm that we have used 
in our mesh generation is an incremental insertion algorithm given by Watson. A brief 
outline of it is stated below : 
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• insert new site Q into existing Delaunay triangulation. 

• find any triangle with circumcircle containing site Q. 

• perform a search to find remaining set of triangles with circumcircle containing site 

Q. 

• construct list of edges associated with deletable triangles. Delete all edges from 
the list that appear more than once. 

• Connect remaining edges to site Q and update Delaunay data structure. 

After obtaining an initial mesh the elements lying in the holes and convex regions axe 
being removed. Finally, based on height functions (based on spacing of the points of the 
perimeter) we insert more elements, where desired, to refine the mesh. 

An automatic unstructured mesh generator uses some sort of boundary data alone 
as an input, and then create the interior nodes and elements automatically based on the 
boundary data. As the mesh generator makes little or no assumption on the structure 
of the domain any structure can be modeled. 

The basic data structure of the mesh consists of X-array and IEN-array. The X- 
array i.e. x(nsd,numnp) which stores the x and y coordinates of the nodal points. Here, 
nsd stands for number of space dimension. It is equal to 2 for our case, numnp is the 
global node number. 

IEN-array, i.e. ien(nen, numel) which defines the connectivity between the global 
nodal points and it contains the information of element numbering with respect to local 
and global node numbering. Here, nen stands for local node numbers, 3 for triangular 
elements and 4 for quadrilateral elements, numel is the global element number. 


3.1 Mesh Design 


Mesh design for flow computation is the most important issue, it must take into account 
the nature of flow. For example, the flow near the body involves boundary layer devel- 
opment, boundary layer separation, formation of shear layer, interaction of shear layers, 
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sharp velocity gradients etc., whereas far away from the body upstream and downstream 
flow changes are small. Around the body we generate structured mesh which is obtained 
by drawing lines perpendicular to the slope of the airfoil surface. This feature is also 
essential from the point of view of turbulence modeling. As boundary layer velocity 
profile does not vary lineaxly, a provision has been made to vary the spacing between 
different layers of the structured mesh in geometric progression. To capture the wake a 
fictitious boundary has been inserted so that the concentration of elements increase in 
the zone near the airfoil. 

Automatic mesh generator description: The grid generator routine needs 
various parameters like the points on the surface of the airfoil, the points on the outer 
domain, the points on a fictitious boundary if present, the first element thickness for the 
structured mesh, the number of elements on the structured mesh and the total thickness 
of structured mesh. For the mesh around multi-element airfoil we need to supply these 
inputs for each element. 

To obtain meshes at certain angle of attack with the free stream velocity, the nodes 
are rotated about certain point, i.e. input parameter. The new coordinates are being 
obtained using the following simple equations. 

x' = xcoscj) — ysin<t > , y' = xsincj) + ycoscj) (3.1.1) 

The mesh generator uses a mesh moving technique such that mesh can be obtained on 
an airfoil at different angles of attack restoring the number of elements and the number 
of nodes. 


3.2 Governing Equations 


Let Cl c R n,d and (0,T) be the spatial and temporal domains respectively, where n s a 
is the number of space dimensions, and T denote the boundary Cl. The spatial and 
temporal coordinates are denoted by x and t. The Navier-Stokes equations governing 
incompressible fluid flow are 


p ( — + u • Vu - f J — V * <r = 0 on Cl for (0, T) 


(3.2.2) 
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V • u = 0 on Cl for (0, T) (3.2.3) 

Here, p, u, f and a are the density, velocity, body force and the stress tensor, respectively. 
The stress tensor is written as the sum of its isotropic and deviatropic parts : 

a — —pi + T , T = 2Mu) , e(u) = |((Vu) + (Vuf) (3.2.4) 

where p and p, are the pressure and coefficient of dynamic viscosity, respectively. Both 
the Dirichlet and Neumann type boundary conditions are accounted for, represented as 

u = g on T g , n • a = h on T h (3.2.5) 

where, T g and are complementary subsets of the boundary T. The initial condition 
on the velocity is specified on Cl : 


u(x, 0) = u o on O 


(3.2.6) 


where, uo is divergence free. The Spalart-Allmaras (one equation turbulence model) for 
eddy viscosity is 


dv 


— + u • Vi/ - c 61 [1 - f t 2] Su - — [V • ((z/ + v) Vi/) + c b2 (Vi/) 2 ] 


+ 


Cwlfw n ftS. 


K 


V 

d J 


— /uAC/ 2 = 0 on Cl for (0,T) 


(3.2.7) 


Here, u is the molecular viscosity and v obeys the above transport equation. So, the 

3 _ 

eddy viscosity [v t ) is defined as v t = vf v i where f vX = , X = ^ > S = S + 

-yjy- , f v2 = l— 1 + *f vl where S is the magnitude of the vorticity, and d is the 

distance to the closest wall. The function f w = 


i+C 


s 6 +c: 


p 


'm3 


] 6 , P = r + c u , 2 (r 6 -r) , r = 
The ft 2 function is given by / t2 = Ct 3 exp (— CuX 2 ) and the trip function / tl is 
as follows : ft 1 = c t igtexp (— c t 2^2 [d 2 + 9 t^t]) where d t is the distance from the field 
point to the trip, which is on the wall, ui t is the wall vorticity at the trip, and A U is 
the difference between the velocity at the field point and that at the trip.Then g t — 
min (o.l, ^ Axt ) where Ax t is the grid spacing along the wall at the trip, whereas the 
others parameters are defined as c b 1 = 0.1355 , & = | , c b2 = 0.622 , k = 0.41 , c^i = 
( 1 +c , 62 j ^ ^2 — 0.3 , c w 3 = 2 , Cyi = 7.1 , Cti = 1 , Ct 2 = 2 , Q3 = 1.2 and Ct 4 = 0.5 
Both the Dirichlet and Neumann type boundary conditions are accounted for, represented 
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as 


v = gi onTg , n - v = h x onY h 


(3.2.8) 


where, T g and F^ are complementary subsets of the boundary T. The initial condition 
on the eddy viscosity is specified on Q, : 

P(x, 0) = D 0 on fl (3.2.9) 


where, v Q is divergence free. 


3.3 Finite Element Formulation 


Consider a finite element discretization of Yl into sub domains f l e ,e — 1, 2, ..., n^, where 
n e i is the number of elements. Based on this discretization, for velocity, pressure and eddy 
viscosity we define the finite element trial function spaces <S£, S g and S P h, and weighting 
function spaces V„, Vp and V P h. These function spaces are selected, by taking the 
Dirichlet boundary conditions into account, as subsets of [H 1A (Q)] nad and H 1/l (f2), where 
H 1/l (Q) is the finite dimensional function space over fl. The stabilized finite element 
formulation of equation 3.2.2 and equation 3.2.3 is written as follows: find u h € <S£ and 
ph 6 gh suc i 1 t h at Vw h e V„, q h G Vp 



+ u h ■ Vu ft - f 


d,Q,+ f e(w /l ) : cr(p h , u h )dYl + f q h V • u h dQ 
Jn Jci 

+ X] / ~ ysuPGpd 1 • Vw 71 + tpspgV q h ) * 


^ ‘ VUk ~ f) “ 


dn e 


n el /• p 

+ Y I 5V -w h pV -u h dQ e = w h • h h dF(3.3.10) 

e=l Jt * Jt i> 


The equation 3.2.7 is written as follows: find v h 6 Sj? h such that Vw ft € V? h and 
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F = c bl [1 - f t2 ] SD h - c b2 (VD h ) 2 + [c wl f w - c $rf t2 ] [^] 2 - f a AU 2 

L v,h -{ ? w +uKv ^- F ) da+ Lv' t Vwl • ((" ■ + ph ) vp ' 1 )} dn 


/n 

n e i 


+±L (tsupgU 11 • Vw A j • 


dv h 

dt 


+ u h • Vir — F 


dQ e 


= [ (3.3.11) 


In the variational formulation, the first three terms of equation 3.3.10, the first two terms 
of equation 3.3.11 and the right hand side of both equations constitute the Galerkin 
formulation of the problem. The first series of element level integrals are the SUPG 
and PSPG stabilization terms added to the variational formulations. In the current 
formulation tpspg is the same as tsupg and is given as 


2Huj 12* 

U h 1 { h 2 


)T 


(3.3.12) 


The second series of element level integrals are added to the formulation for numerical 
stability at high Reynolds numbers. This is a least squares terms based on the continuity 
equation. The coefficient 5 is defined as 


<5 = 



(3.3.13) 


where, 


2 = 


(^) Re u < 3 
1 Re u > 3 


(3.3.14) 


and Re u is the cell Reynolds number. Both stabilization terms are weighted residuals, 
and therefore maintain the consistency of the formulation, h is the element length and 
various definitions of h is available. The one which results in the least sensitivity of the 
computed flow to the element aspect ratio has been used for computations in the present 
work. According to this definition, the element length is equal to the minimum edge 
length of a triangular (3 noded) element. 


The time discretization of variational formulations given by the equations 3.3.10 
Sz 3.3.11 is done via the generalized trapezoidal rule (Crank-Nicholson). For unsteady 
computations, we employ a second order accurate in time procedure. Equal in order basis 
functions for velocity and pressure(the P1P1 element) are used and a 3 point quadrature 
is employed for numerical integration. The non-linear equation systems resulting from 
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the finite element discretization of the flow equations are solved using the Generalized 
Minimal RESidual(GMRES) technique. The pre-conditioning used in conjunction with 
the GMRES method is nodal block diagonal. Matrix free GMRES algorithm is used 
which reduces the memory requirements. The key advantage of the matrix free GMRES 
algorithm is direct computation of the result of matrix vector products in the GMRES 
algorithm which avoids the explicit formation of the element level matrices. 



Chapter 4 


Results and Discussions 


This chapter presents the results of 2-D simulations of turbulent flows using Spalart- 
Allmaras turbulence model on three geometries of increasing complexity. Initially the cal- 
culations have been done for flat plate, then NACA 0012 airfoil and finally for McDonnell 
Douglas 30P-30N multi-element airfoil. The incompressible Reynolds Averaged Navier 
Stokes equations, in conjunction with the Spalart-Allmaras turbulence model for turbu- 
lence closure, have been solved using stabilized finite element formulations. To stabilize 
the computations against spurious numerical oscillations in advection dominated flows 
and to enable the use of equal order interpolation velocity-pressure elements, streamline- 
upwind/Petrov-Galerkin(SUPG) and pressure-stabilizing/Petrov-Galerkin(PSPG) stabi- 
lization techniques have also been employed. In this technique, stabilizing terms are 
added to the basic Galerkin formulation (Chapter 3). Equal-in-order linear basis func- 
tions for velocity and pressure have been used and a 3 point quadrature has been em- 
ployed for numerical integration. The time integration of the equations has been carried 
out via the time accurate version of the generalized trapezoidal rule(Crank-Nicholson). 
The nonlinear implicit equation system resulting from the finite element discretization of 
the flow equations has been solved using the Generalized Minimal RESidual (GMRES) 
technique in conjunction with diagonal preconditioners. 

The finite element meshes used in various computations consist of a structured mesh 
close to the body and an unstructured mesh generated via Delaunay’s triangulation, in 
the rest of the domain. This type of grid has the ability of handling fairly complex 
geometries while still providing the desired resolution close to the body to effectively 
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capture the boundary layer and shear layer, especially in the context of unsteady flows. 
In addition, the unstructured grid around the body allows for efficient implementation 
of the turbulence model. 

The non-dimensional quantities encountered in discussing the results are : 


• Pressure coefficient. 


• Drag coefficient. 


• Lift coefficient. 


• Moment coefficient. 


C p = 


(g ~ Poo) 
jjPooll^OolP 



where, A is the reference area, l is the reference length, D is the drag force, L is 
the lift force and M is the pitching moment. 


4.1 Flat Plate 

4.1.1 Computational Domain & Boundary conditions 

The flat plate resides in a rectangular computational domain of dimension 15x5 units. 
The inflow boundary is located at a distance of 5 unit from the leading edge of the plate 
as shown in Figure 4.1. On the plate surface, the no-slip condition for the velocity and 
zero value for eddy viscosity have been specified while free stream values have been 
assigned for the velocity at the inflow boundary, whereas a very low value of eddy 
viscosity (P = 10~ 10 ) has been specified at the same boundary. At the downstream 
boundary, a Neumann type boundary condition for both velocity and eddy viscosity 
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has been specified. The Reynolds number is calculated based on unit length of plate, 
free-stream speed and viscosity of the fluid. 


4.1.2 Turbulent calculation for different values of Re 


Results have been presented in this section for different values of Re using two different 
Meshes (one coarse & one fine), as shown in Table 4.1. 


Re 

Meshl 

nodes elements 

Mesh2 

nodes elements 

10 s 

5151 

5000 

15251 

15000 

10 6 

5151 

5000 

15251 

15000 


Table 4.1: Number of nodes and elements in different meshes used for turbulent calcu- 
lation of flow past a flat plate. 


The turbulent computations have been carried out in two-dimensions by specifying 
the transition location on the plate surface for Re = 10 5 & 10 6 . The variation of 
skin-friction coefficient^ f) , turbulence index(it ) along the length of the plate have been 
shown. The computed results have also been compared with the theoretical values of 
(C f ), calculated from correlation function given by F. M. White [38], defined as: 

f laminar = > 

C f = \ turbulent - - 

1 turouiem - Jn 2 (0 . 06Rei ) 

where x is the distance along the length of the plate. 

Figure 4.2 & Figure 4.3 shows the variation of Cj along the length of the plate 
for Re = 10 5 & 10 6 , respectively. It is evident that the computed results are in good 
agreement with the theoretical values. Also, the desired trip location, where the flow 
becomes turbulent, can be controlled by specifying the trip point on the plate surface. 
One interesting phenomenon observed here is that, as we keep on increasing the Re, the 
trip location starts moving towards the leading edge. So in that case if the trip point 
is specified towards the trailing edge, no effect of the specified point is sensed by the 
flow field and it trips on its own at an upstream location. Same features have been 
noticed even after the spatial resolution increased to a sufficient level as used in Mesh2. 
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Figure 4.4 shows the case where the flow trips on its own at an upstream location 
due to increase of Re from 10 5 to IQ 6 even after the desired trip point is specified. 

Figure 4.5 shows the variation of turbulence index along the length of the plate. It 
is seen from the figure that at the specified trip location, the turbulence index achieves 
an approximate value of 1 and sustains this value, thereafter. In Figure 4.6, the turbu- 
lent boundary layer velocity profile is shown, which is also in good agreement with the 
theoretical value. 
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Figure 4.1: Flat plate configuration. 
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Figure 4.2: Variation of C f along the length of the plate for Re = 10 5 . 
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4.2 NACA 0012 airfoil 

4.2.1 Computational Domain &; Boundary conditions 

The NACA 0012 airfoil lies in a rectangular computational domain whose upstream and 
downstream boundaries are located at 5 and 11 chord lengths from the leading edge, 
respectively. The upper and lower boundaries are placed at 5 chord lengths, each, from 
the leading edge. On airfoil surface, no-slip boundary condition for velocity and zero 
value for eddy viscosity have been specified while free-stream values have been assigned 
for both velocity and eddy viscosity at the upstream boundary. At the downstream 
boundary, a Neumann type boundary conditions for both velocity and eddy viscosity 
has been specified that correspond to zero viscous stress vector and zero eddy viscosity, 
respectively. At the upper and lower surface boundaries, the component of velocity 
normal to the component of stress vector along these boundaries have been prescribed 
zero value. The Reynolds number based on the chord length of the airfoil, free-stream 
velocity and viscosity of the fluid is 10 6 . 


4.2.2 Finite element mesh and Results 


The computed results have been reported in this section, carried out using typical finite 
element mesh as shown in Figure 4.7. The unstructured mesh provides flexibility to 
handle complex geometries. The structured mesh around the airfoil provides effective 
control on the grid to resolve the boundary layer. The meshes used for the turbulent 
calculation at different angle of attack have been listed in the Table 4.2. 


a 

Mesh 

nodes elements 

5° 

18925 

37604 

10° 

18994 

37742 

11° 

18994 

37742 

15° 

19071 

37896 


Table 4.2: Number of nodes and elements in NACA mesh at different angle of attack. 


The number of nodes and elements for a = 10° k a = 11° are same as the mesh 
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for a — 11° is obtained by deforming mesh for = 10° using mesh-moving scheme. 

The initial condition for the computation for = 1 1° is the fully developed solution 
for a = 10°. Figure 4.8 shows the pressure distribution over airfoil surface at different 
angle of attacks. Also, It shows the comparison with Baldwin-Lomax solution without 
any significant deviation in the result. At various angle of attacks, the velocity and 
eddy viscosity plots are shown in Figure 4.9. It is noticed that the flow is attached upto 
a = 11° and starts separating at a = 15° towards the trailing edge of the airfoil where 
the boundary layer thickness is thicker than the leading edge whereas the thickness on 
the lower surface is much much thinner than that on upper surface due to the difference 
in the pressure gradients. The variation of Ci,C d ,C m are shown in Figure 4.10 which 
also represents the same phenomenon that the flows started separating at a = 15° where 
the C d value increases significantly. In Figure 4.11, the variation of turbulence index over 
airfoil surfaces is shown for different angle of attacks. The turbulence index variation 
also estimates the same flow behaviours that the flow starts separating at a = 15°. The 
chord wise locations for velocity profiles, shown in Figure 4.12, are listed in the Table 4.3. 


Station No 

x/c 

surface 

1 

1.00 

lower 

2 

0.240 

lower 

3 

0.052 

lower 

4 

0.008 

upper 

5 

0.207 

upper 

6 

1.00 

upper 


Table 4.3: NACA 0012 airfoil: Location of stations from leading edge for showing velocity 
profiles. 


Figure 4.13, Figure 4.14 & Figure 4.15 shows the nature of the velocity profile at 
different stations over airfoil surfaces for various angle of attacks. The computed results 
have also been compared with the results, obtained using Baldwin-Lomax turbulence 
model. It shows that the computed results are in good agreement with Baldwin-Lomax 
solutions. 
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Figure 4.7: Finite element mesh for NACA 0012 airfoil at a = 10°. 
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x/c 



x/c 



Figure 4.8: C p distribution over NACA 0012 airfoil surface at different angle of attacks 
for Re = 10 6 . 






4.2 NACA 0012 airfoil 


39 



Figure 4.9: Turbulent flow past a NACA 0012 airfoil at Re - 10 6 : magnitude of veloc- 
ity(left) and eddy viscosity(right) for various angle of attacks. 
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station 4(x/c = 0.008) 


station 5(x/c = 0.207) 


station 6(x/c =1.00) 




station 3(x/c = 0.052) 


station l(x/c = 1.00) 


station 2(x/c = 0.240) 


Figure 4.12: NACA 0012 airfoil: Location of stations from the leading edge where velocity- 
profile, over airfoil surface is shown. 
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Figure 4.14: NACA 0012 airfoil: Velocity profiles at different stations over airfoil surface 
for a — 10°. 
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4.3 McDonnell Douglas Multi-Element Airfoil 

4.3.1 Airfoil configuration 

The high lift model investigated is a McDonnell Douglas aerospace (MD A) 2D 11.55% 
thick super critical, single flap, three element airfoil as shown in Figure 4.16(a). The slat 
chord is 14.48% and the flap chord is 30% of the stowed airfoil chord. The airfoil was 
configured in a typical approach/landing configuration with slat and flap deflections of 
30°. This configuration is denoted as 30P-30N(slat deflection 30°, slat rigging positive, 
flap deflection 30°, flap rigging negative). Figure 4.16(b) defines the nomenclature for 
gap and overhang. The slat and flap setting for 30P-30N configuration has been shown 
in Table 4.4. 


Nomenclature 
(Flap & Slat) 

MD Geometry 
(30P-30N) 

Slat Deflection((5 s ) 

-30° 

Slat Gap, %c 

2.95 

Slat Overhang, %c 

-2.5 

Flap Deflection (5f) 

30° 

Flap Gap, %c 

1.27 

Flap Overhang, %c 

0.25 


Table 4.4: Slat and Flap setting of Multi-Element airfoil 


4.3.2 Computational Domain and Boundary Conditions 

The airfoil rests in a rectangular domain. The grid extends 7c upstream and 11c down- 
stream of the model to minimize inflow/outflow influence on the solution. The details of 
the meshes used for computations for various angle of attack have been summarized in 
Table 4.5. A close-up of the grid is shown in Figures 4.17 and 4.18. The structured mesh 
thickness over each body varies from 0.1c to 0.007c. The meshes used in computation 
contains 457 points on main element, 198 points on slat and 200 points on flap. In order 
to control the unstructured mesh density in wake region, three fictitious boundaries are 
inserted to resolve the wake effectively. The spacing of the first point normal to the wall 
is is 1 x 10 -6 normalized to the chord length of the airfoil in the undeflected position. 
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This spacing yields a y + of less than 1.0 for the point next to the wall over main element. 

On airfoil surface, no slip boundary condition for velocity and zero value for eddy 
viscosity have been specified while free stream values have been assigned for both velocity 
and eddy viscosity at the upstream boundary. At the downstream boundary, a Neuma n n 
type boundary condition for both velocity and eddy viscosity has been specified. At the 
upper and lower surface boundaries, the component of velocity normal to the component 
of stress vector along these boundaries have been prescribed zero value. 


a 

Meshl 

nodes elements 

Mesh2 

nodes elements 

8° 

46916 

92927 

84320 

167735 

~w~ 

48760 

96615 

- 

- 

wl 

48760 

96615 


- 

HF 1 

48368 

95831 

- 

- 

21° 

48236 

95567 

83052 

165199 

27° 

47669 

94433 

- 

- 


Table 4.5: Number of nodes and elements in Multi-Element mesh at different angle of 

attack. 

The number of nodes and elements for a = 15° & a = 16° are the same since the 
mesh for a = 16° is obtained by deforming the mesh for = 15° using a mesh-moving ' 
scheme. 


4.3.3 Results 


This section contains the comparison of computed and experimental results for various 
angle of attack at Re = 9 x 10 6 . For both the meshes the first point next to the wall 
is at a distance of 1.0 x 10" 6 , which gives an average y + of 0.5 on the upper surface of 
main element. Mesh2 has also been used for computation at a = 8° & a = 21° to study 

mesh convergence. 


4. 3.3.1 Pressure distribution 
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a = 8 ° 

Figure 4.19 shows the comparison of computed and experimental pressure distribu- 
tion over slat, main element and flap for a = 8°. It has been observed that the computed 
results for meshl predict higher values whereas mesh.2 predict lower values than the ex- 
perimental values on the upper of the slat. However, on the lower surface of the slat, 
the computed results are in good agreement with the experimental results. 

1 he main element pressure distribution shows a large suction peak near the leading 
edge followed by pressure recovery region. The pressure is not recovered to free stream 
level since the trailing edge is adjacent to the flap suction peak. The lower surface is 
characterized by high pressure (nearly stagnant flow C p = 1). The free stream pressure in 
the present case has been assigned to zero for the computation of C p distribution. On the 
upper surface, the computed results for meshl are in good agreement with experimental 
results, whereas the results for mesh2 predict slightly lower values than the experimental 
values. 

The flap upper surface pressure distribution shows higher pressure recovery for 
meshl results than experiment, indicating higher separation near the trailing edge. How- 
ever, the results obtained using SA model for mesh2 are in good agreement with the 
experimental results, while BL results for mesh2 show lesser pressure recovery than ex- 
periment, indicating lesser separation near trailing edge. 


a = 15° 


Pressure distribution on the slat, main element and flap has been shown in Fig- 
ure 4.20 for a = 15°. No experimental data is available for this case, so the results have 
been compared for two different turbulence models. On the upper surface, BL shows 
higher suction peak near the leading edge of slat and flap, indicating higher pressure 
recovery. While both the models predict same results for main element. 

a =16° 


Figure 4.21 shows pressure distribution on the slat, main element and flap for 
a = 16®. The computed results are in good agreement with the experimental results, 
except on the upper surface of the flap. On the flap upper surface, the SA model shows 
higher pressure recovery than experiment. Also it predicts larger suction peak than 
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experiment near the leading edge of the main element, 
a = 19°, Re — 7x 10® 


Pressure distribution on the slat, main element and flap for a = 19° at Re = 7x 10 6 
has been shown in Figure 4.22. Experimental data is not available for this case, so the 
results have been reported for two different turbulence models. On the upper surface, 
BL predicts higher suction peak near the leading edge of slat and flap, indicating higher 
pressure recovery. While both the models show same results for main element. 

a = 19°, Re - 9 x 10* 

Figure 4.23 shows pressure distribution on the slat, main element and flap for 
Q! =: 19° ■h’- e = 9 x 10 6 . The suction peaks are larger in this case than that at lower 

angle of attack. Again there is not much difference in pressure distribution for computed 
and experimental results, except on the upper surface of the flap. Near the leading edge 
of main element, SA model predicts larger suction peak than experiment. The computed 
pressure distribution of SA model is more close to the experiment than BL model in all 
the three elements. 


a = 21 ° 


Figure 4.24 shows pressure distribution on the slat, main element and flap for 
a = 21°. It has been observed that the computed results predict higher suction peak 
and also show significant differences from the experimental results on the upper and lower 
surfaces of the slat. However, on the main element the computed pressure distribution 
matches well with the experiment. While showing the C p distribution over flap surface, 
SA model shows good agreement with the experiment, whereas BL results for both 
meshes predict higher values than experiment, followed by higher pressure recovery, 
indicating higher separation near trailing edge. 


a = 27° 


Pressure distribution on the slat, main element and flap has been shown in Fig- 
ure 4.25 for a — 27°. Experimental data is not available for this case also, so the results 
have been compared for two different turbulence models. The suction peaks are having 
large values in this case than that at lower angle of attack. BL model shows larger sue- 
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tion peak near the leading edge of slat and flap than SA model, whereas in case of main 
element it shows same value of suction peak with lesser pressure recovery on the upper 
surface. The flap upper surface pressure distribution shows higher pressure recovery for 
BL model than SA model, indicating higher separation near the tr ai ling edge. 

In general, the agreement between the C p data of the experiment and all two of 
the turbulence models is quite good. The biggest discrepancies occur for the lower 
angle of attack, particularly on the slat and the flap. The greater the region of flow 
separation on the flap, the lower the lift on the flap and slat. The BL model shows 
higher flow separation on the flap: a flattening of the C p distribution near the trailing 
edge. It has been observed that the SA model results in C p distributions that are closer 
to experimental observations. 


4.3.3. 2 Velocity profiles 


Velocity profiles have been presented at six stations for a = 19° and four stations for 
rest of the angle of attacks. The chord wise locations for velocity profiles, shown in 

Figure 4.26, are listed in the Table 4.6. 


Station No 

x/c 

1 

0.1075 

2 

0.45 

3 

0.85 

4 

0.89817 

5 

1.0321 

6 

1.1125 


Table 4.6: Multi-element airfoil: Location of stations from leading edge for showing 

velocity profiles. 


a = 8° 

Figure 4.27 k 4.28 show the contours of the velocity magnitude, pressure fields, 
eddy viscosity for meshl and mesh2, respectively. Velocity profiles at different stations 
have been shown in Figure 4.29. According to the experimental data, the main element 
boundary layer edge is just below n/c = 0.01 (Fig. 4.29, station 1). The slat wake passes 
just above the boundary layer and shows very little velocity deficit. The result obtained 
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using BL model for meshl is in good agreement with the experimental result, whereas 
the BL result for mesh.2 shows higher momemtum deficit in the slat wake and thicker 
boundary layei than the experiment. The SA result for meshl is very close to the 
experimental result but shows slightly wider slat wake. The SA result for mesh2 does 
not show good match with experimental results. 

Velocity profile on the flap near the leading edge has been presented in Fig- 
ure 4. 29 (station 4). The profile in this region usually consists of four levels. Going from 
the surface outward these four levels are: (1) flap boundary layer, (2) slot flow through 
the flap gap, (3) main element wake, and (4) slat wake. Slat boundary layer is extremely 
thin at this location. Experimentally, the slot flow extends from 0.001 < n/c < 0.01 and 
is ramp shaped indicating strong viscous effects on the slot flow coming through the flap 
gap. These viscous effects are quite possibly generated by the separated and recirculating 
flow in the flap cove. The main element wake extends from n/c = 0.01 to n/c = 0.025. 
Both the models for meshl show good agreement with the experimental result while BL 
model for mesh2 under predicts the edge velocity in flap boundary layer and also predicts 
too deeper slat and main element wake and no where close to the experimental result. 
Also, the SA model for mesh2 under predicts the velocity magnitude above the main 
element wake and not in good agreement with the experiment. Experimental result does 
not show any slat wake deficit whereas computed results show wide slat wake. 

The velocity magnitude obtained using SA model for both meshes and BL model 
for meshl is almost close to the experimental result while BL model for mesh2 predicts 
too deeper main element and slat wake(Fig. 4.29, station 5). The velocity magnitude of 
flap boundary layer near trailing edge(Fig. 4.29, station 6) is completely dispersed. 


a = 15° 


The velocity magnitude, pressure fields, eddy viscosity contours have been shown 
in Figure 4.30. Velocity profiles on main element for all three stations have been shown 
in Figure 4.31 and on flap in Figure 4.32 for two different turbulence models. There are 

no experimental results available for this angle of attack. 

a — 16° 

Figure 4,33 shows the velocity magnitude, pressure fields, eddy viscosity contours 
for 16° angle of attack. Velocity profiles at various stations have been shown m Fig- 
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ure 4.34. At the leading edge of the Deputation 4), the velocity profile shows good 
agreement with experiment in predicting main element walce. However, on main .w™* 
it shows higher velocity deficit in the slat waie, but it predicts overall good result with 
experiment. While considerable scattering has been observed in the prediction of flap 
velocity profiles, especially at the station 5(mid-chord position) and station 6(trailing 
edge) of the flap. 

a = 19°, Re = 7 x 10® 

Velocity magnitude, pressure fields, eddy viscosity contours have been shown in 
Figure 4.35. Velocity profiles on main element for all three stations have been shown in 
Figure 4.36 and on flap in Figure 4.37 for two different turbulence models. There are no 
experimental results available for this angle of attack. 

a = 19°, Re = 9 x 10® 

Figure 4.38 shows the velocity magnitude, pressure fields, eddy viscosity contours 
for a = 19°, Re = 9 x 10°. Velocity profiles on main element for all three stations have 
been shown in Figure 4.39. In this case slat wake has been predicted nicely well by both 
the turbulence models. The results obtained by SA model are in good agreement with 
the experimental results, while BL model over predicts the results. Velocity profiles on 
the flap have been presented in Figure 4.40. At the leading edge of the flap(station 4), the 
velocity profile shows good agreement with experiment in predicting main element wake. 
While SA model shows higher velocity deficit in the slat wake, but it predicts overall 
good result with experiment. There is considerably more scattering in the prediction of 
flap velocity profiles than those of main element, especially at the station 5(mid-chord 
position) and station 6(trailing edge) of the flap. 


a = 21 ° 


The velocity magnitude, pressure fields, eddy viscosity contours have been pre- 
sented in Figure 4.41, 4.42 for meshl k mesh2, respectively. Velocity profiles at different 
stations have been presented in Figure 4.43. In this case slat wake has been predicted 
nicely by SA turbulence model. The SA model results show overall good agreement 
with the experimental results. On main element (station 1), the SA result for meshl 
quite close to the experimental result but shows slightly higher velocity deficit in the slat 
wake, while the SA result for mesh2 predicts much more higher velocity deficit in the 
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slat wake. At the leading edge of the flap(station 4), the velocity profile computed by 
SA model shows good agreement with experiment in predicting main element wake as 
well as exhibits overall good agreement. While significant scattering has been observed 
in the piediction of flap velocity profiles, especially at the station 5(mid-chord position) 
and station 6 (trailing edge) of the flap. 


a = 27° 


Figure 4.44 presents the velocity magnitude, pressure fields, eddy viscosity contours 
for ex = 27°. Velocity profiles on main element for all three stations have been shown in 
Figure 4.45 and on flap in Figure 4.46 for two different turbulence models. There are no 
experimental results available for this angle of attack. 

In general, the agreement between the experimental and computational velocity 
profiles is fairly good. One of the biggest differences between the computations has been 
observed at a = 8°: the velocity defect of slat wake is very small(at the station 1) in 
the experimental data, whereas all the computations predict much larger defect. This is 
bit larger in mesh2 results, and less so on the meshl calculations. So, it is difficult to 
conclude which model is giving the best results overall, since none of the models agree 
completely with the experimental results, and differences between different computations 
are small. 

It is well known that stall occurs near a = 21°. But it is clear from Figure 4.44 
that the flow over main element and flap is still fully attached and as a result of that 
both the models are not capable of predicting stall even at a = 27° . 


4.3.3. 3 Force and Pitching moment 

Computed lift force(at Re = 9 x 10 6 ) has been compared with experimental results and 
presented in Figure 4.47. The experimental lift curve shows relatively linear behaviour 
up to a. — 12 °, followed by a gradual rounding, whereas Ci max occurs near a — 21 and 
mildly stalls thereafter. Computed results also show similar trends but predict higher 
Cimax than the experimental result. The lift curve obtained by SA model shows linear 
behaviour upto a = 15°, and thereafter looses its linearity. The flap lift compares very 
well with the experiment. While the slat lift is slightly higher than experiment and the 
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main element lift compares sufficiently good before ce = 19° but the experimental curve 
shows a drop in lift at a = 21°. There seems to be a constant increment in lift in the 
computed results. The lift values are all quite close up to 16°, as anticipated by the C p 
results. None of the models agree with experimental value of maximum lift. This is most 
likely because the experiment does start to undergo some 3D effects become sig nifi cant 
at high values of lift. 

The drag force has been shown in Figure 4.48. It has been observed that Cd is 
not well predicted by any models but SA model’s prediction is better than BL model. 
In the present work the drag has been computed by integrating the pressure and skin 
friction forces on the surface. This method has been shown [36] to be less accurate than 
other methods(such as wake integration), and may be the cause of large errors in the 
computed C,i . Pitching moment has been shown in Figure 4.49. The pitching moment 
coefficient predictions are good for SA model at the lower angle of attack, although they 
differ from BL model. The trend of increasing slope in C m between a = 12° and 16° is 
not predicted by any model. Comparison shows that SA model over predicts the result 
while BL model under predicts. 
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Figure 4.16: McDonnell Douglas Multi-Element Airfoil: Model geometry & Nomencla 

ture (as given in the paper [17]). 



4.3 McDonnell Douglas Multi-Element Airfoil 


56 





Figure 4.17: Typical Finite Element Mesh shown for a — 15°, consists of 48760 nodes 
and 96615 elements. 





Figure 4.18: Finite Element Mesh for a = 15° showing close-up of main-element, slat 
and flap. 
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Figure 4.19: Pressure distribution over multi-element surfaces for a = 8°. 
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Figure 4.20: Pressure distribution over multi-element surfaces for 
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Figure 4.22: Pressure distribution over 


multi-element surfaces for a = 19°, Re = 7x 10 6 . 




4.3 McDonnell Douglas Multi-Element Airfoil 


62 


slat 



main 


o. 

o 


rrprr i t-pnt t i ; ni - rp rrrjTT i i | rr r Tj 

Spalart-Allmaras 

Baldwin-Lomax 

exp o 



t — © — 0- « © 0 0 © — tt&tKrOWrv 

ElJU ,i t..l i. iiJ J 


0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 
x/c 


flap 



Figure 4.23: Pressure distribution over multi-element surfaces for a = 19°, Re - 9 x 10 6 






4.3 McDonnell Douglas Multi-Element Airfoil 


63 



-0.10 -0.08 -0.06 -0.04 -0.02 0.00 0.02 


x/c 


main 



0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 


x/c 



Figure 4.24: Pressure distribution over 


multi-element surfaces for a = 21°. 
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Figure 4.25: Pressure distribution over multi-element surfaces for a 27° 
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Figure 4.26: Multi-Element Airfoil: Location of stations from leading edge, for showing 

velocity profiles. 
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Figure 4.28: Multi-element airfoil: Contours of velocity magnitude, pressure fields, eddy 
viscosity for rnesh2 at a = 8°, Re = 9 x 10 6 . 
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Figure 4.29: Multi-element airfoil: velocity distribution at a = 8°, Re - 9 x 10 6 
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Figure 4.31: Multi-element airfoil: velocity distribution on mam element at a 15 

Re = 9 x 10 6 . 
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Figure 4.32: Multi-element airfoil: velocity distribution on flap at a - 15°, He - 9 x 10 6 
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Figure 4.35: Multi-element airfoil: Contours of velocity magnitude, pressure fields, eddy 

viscosity at a = 19°, Re = 7 x 10 6 . 
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Figure 4.36: Multi-element airfoil: velocity distribution on main element at a - 19° 

Re = 7 x 10 6 . 
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Figure 4.37: 
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Multi-element airfoil: velocity distribution on flap at a - 19°, Be - 7x 10 6 
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Figure 4.38: Multi-element airfoil: Contours of velocity magnitude, pressure fields, eddy 
viscosity at ot — 19°, Re = 9 x 10 6 . 
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Figure 4.39: Multi-element airfoil: velocity distribution on main element at 
Re — 9 x 10 6 . 


= 19 °, 






4.3 McDonnell Douglas Multi-Element Airfoil 


79 



0.0 0.5 1.0 1.5 2.0 


station 5 



station 6 



Figure 4.40: Multi-element airfoil: velocity distribution on flap at a = 19°, Re- 9 x 10 6 . 
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Figure 4.41: Multi-element airfoil: Contours^of velocity magnitude, pressure fields, eddy 
viscosity for meshl at a = 21°, Re = 9 x 10 . 
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Figure 4.42: Multi-element airfoil: Contours of velocity magnitude, pressure fields, eddy 
viscosity for mesh.2 at a — 21°, Re = 9 x 10 . 
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Figure 4.43: Multi-element airfoil: velocity distribution at a - 21°, JRe 9 x 10 
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Figure 4.44: Multi-element airfoil: Contours of velocity magnitude, pressure fields, eddy 
viscosity at a = 27°, Re = 9 x 10 6 . 
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Figure 4.45: Multi-element airfoil: velocity distribution on main element at a = 27°. 
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Figure 4.46: Multi-element airfoil: velocity distribution on flap at a = 27°. 
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Figure 4.48: Multi-element airfoil: C d Vs angle of attack at Re = 9 x 10 6 
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Figure 4.49: Multi-element airfoil: C m Vs angle of attack at Re - 9 x 10 6 . 
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4.4 Convergence of GMRES solver 

In the present study, GMRES technique in conjunction with diagonal preconditioners 
has been used to solve the equation system resulting form finite element formulations. It 
was observed that for high aspect ratio elements, even with a fully implicit method, the 
program has a problem in convergence with larger time steps. This was also investigated. 

Laminar Flow 

The solver was working properly in the context of laminar flow. The convergence 
rate was also found impressive. 

Turbulent Flow 

In case of turbulent flow, problem aroused in the solver with poor convergence rate. 
So, it is evident that the problem was caused due to the turbulence model. Therefore, the 
term wise investigation was carried out to improve the convergence rate. Equation 2.6.19 
refers the differential equation for full turbulence model. As, the whole computation has 
been done for fully turbulent cases, the trip terms(/ tl & fa) in equation 2.6.19 were 
assigned to 0. So, the resultant equation which was effectively solved, is equation 2.6.17. 
The following terminologies were used during convergence study. 


DD 

Dt 

Terml 

CbiSP 

Term2 

1 

<r' 

V. ((i/ + v ) Vv) + q ,2 (Vi;) 2 


Term3 

C'wlfw 

V 

d 

i 

Term4 

Change of f w w.r.t. v 

Term5 


Combinations of these terms were includes/excluded to form the LHS of the ma- 
trix resulting from the finite element discretization. During convergence study, always 
changes were made in both preconditoner and the iterative subroutine. The results have 
been tabulated below: 
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Terml 

Term2 

term3 

Term4 

Term5 

Remarks 

On 

Off 

Off 

Off 

Off 

Good convergence 

On 

On 

Off 

Off 

Off 

Good convergence 

On 

On 

On 

Off 

Off 

Poor convergence 

On 

On 

Off 

On 

Off 

Good convergence 

On 

On 

Off 

On 

On 

Poor convergence 


After the above said investigation was carried out, it was found that Term3 & 
Term5 were causing problem. Further analysis led to the approximation of the parameter 
{v + v) VP in the preconditioner and finally this parameter was approximated to (v)W 
to improve the convergence rate, whereas Term5 still remained as a problematic term. 
As a result of that Term5 was totally turned off from both preconditioner & the iterative 
subroutine in the rest of the computations for getting improved convergence rate. 

In addition to these changes, several more things like different values for element 
length (fi) & stabilization term(r), addition of least square terms, larger Krylov spaces, 
were incorporated to improve the convergence rate. But the implementation of all these 
changes had a very little effect in convergence rate, so, these changes were no more used 
in the computations. 

After stabilization, the modified solver was working fine for NACA 0012 airfoil but 
once again started causing problem while implementing for multi-element airfoil with 
larger time steps. To address this problem, computations for the multi-element airfoil 
were carried out in block-iteration manner. First the eddy viscosity equation is solved 
and then velocity & pressure equations. The process is repeated for various iteration 
within a time step. 




Chapter 5 


Conclusions 


2-D simulations of turbulent flows using Spalart-Allmaras turbulence model on three 
geometries of increasing complexity, have been carried out by solving incompressible 
Navier-Stokes equation, in the context of finite element formulation, to understand the 
flow physics as well as the numerical aspects. Initially the computations were done on 
flat plate, then NACA 0012 airfoil and finally for McDonnell Douglas 30P-30N multi- 
element airfoil. Flat plate results showed very good agreement with the theoretical 
results with successful control over the trip locations. The computed results showed 
consistency in the context of increased spatial resolution also. In case of NACA 0012 
airfoil, the obtained results have also been found quite impressive. With the discussed 
results, it can be concluded that this turbulence model is able to capture the flat plate 
flow phenomenon and flow over NACA 0012 airfoil successfully. 

The flow around multi-element airfoil is quite complex involving confluent shear 
layers. The flow on each element should be resolved properly since the shear layer of 
an element affects the flow on subsequent elements. The results have been compared 
between two different turbulence models and experiments. The surface C p distribution 
shows very little difference in the flow solutions of each of the models. The velocity 
profiles exhibit some differences in the computations, particularly in the spreading rates 
of wake regions of the flow. The SA model showed almost close match with the experi- 
ment in the wakes, while the BL model consistently showed the least amount of velocity 
defect in wakes. The most significant difference observed between the computational 
results and the experimental data, is in the slat wake at the lower angle of attack. The 
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grid resolution study showed that this difference did not decrease with increased grid 
resolution. Prom the presented results it is not well understood that the discrepancy in 
the slat wake is a problem with the turbulence models, or with the grid resolution, or 
due to some other factor. The lift coefficient has been over predicted in computations. 
The computed drag is not also very much close to the experimental values. Similarly the 
pitching moment has also been over predicted. Both the models are unable to predict 
stall behaviour even at higher angle of attack, whereas experimental results show stall 
near a = 21°. Given the complexity of the multi-element airfoil geometry and the sur- 
rounding flow field, it would be wisest to proceed with additional treatment of transition 
terms in the SA model which may be helpful in predicting the stall properly. 



Chapter 6 


Scope for Future Work 


• The effect of transition location is explored in Spalart-Allmaras model. But in the 
present study, the code does not employ transition modeling. As a result of that 
the fully turbulent solutions have been accomplished. So, the transition modeling 

can be explored. 

• Hysteresis effects can be studied over multi-element airfoil. 

9 dynamics of deployment of slat and flap can be explored. 

• The fully implicit program has a problem while computing on multi-element ge- 
ometry. So, the GMRES solver can be further analyzed for calculating suitable 

preconditioner. 
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